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FOREWORD 


Recent progress in instrumentation for data acquisition and computing methods suggests 
that satellite geodesy must concern itself with even such minute effects as the gravitational 
action that the tidal dislocations of the earth’s masses exert upon an artificial satellite. A 
number of computer algorithms suited to this purpose were developed at the Naval Surface 
Weapons Center (NSWC) during the past few years. The relevant physics background as well 
as sketches of the computer routines were published in a sequence of four technical reports 
in 1976 and 1977. Since then, the algorithms were worked out in detail, coded for the 
computer, and subjected to the customary numerical checkout. Also, an alternate procedure 

was developed for the ocean tides, which increased the number of tide, routines from the 

original four to five. The necessary derivations for the additional routine were never 
documented. They are included as Appendix A. The present report is intended to serve as 
a manual of the five algorithms, each being presented in a form ready for immediate 
computer implementation, complete with checkout data. 

The work documented here was done in the Astronautics and Geodesy Division and 
was funded as part of the development of the TERRA computer program which is used to 

determine the earth’s gravity field from terrestrial and satellite observations. Henry E. 

Castro of the Computer Program Division coded those algorithms that involve the 
atmospheric tides and the two versions of the ocean tide. R. Gordon Barker, also of the 
Computer Program Division, coded the solid earth tide. Additionally, these two gentlemen 
invested a great effort while making the required checkout runs on the computer. Raymond 
B. Manrique of the Astronautics and Geodesy Division and the author jointly outlined the 
computer algorithms, with the author bearing the overall responsibility for the finished 
formulations and the sole responsibility for the numerical checkout calculations. 

This report was reviewed by Richard J. Anderle, Head, Astronautics and Geodesy 
Division. 


Released by: 

R. A. NIEMANN, Head 
Strategic Systems Department 
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INTRODUCTION 


Several years ago, while compiling a force model on which to base the trajectory 
integration in the TERRA 1 and CELEST 2 computer programs for satellite geodesy and 
satellite ephemeris computation, the author faced the nearly complete lack of published 
algorithms that would in a suitable form express the gravitational action by which the tidal 
dislocation of the earth’s masses perturbs satellite motion. At the same time, it was realized i 

that the earth tides would have to be taken into account if TERRA and CELEST were to 
substantially exceed the accuracy requirements of their forerunners, GEO 3 and ASTRO. 4 i 

An effort was consequently made to derive the necessary disturbing potentials and ^ 

perturbing accelerations. The result consisted of a number of mathematical models, each 
designed to be realistic and describing one type of perturbation corresponding to one of 
the four earth tides. The details of these models were documented in four technical 
papers. 5,6,7,8 Simultaneously, the necessary computer routines were developed and 

submitted for coding. 

The first tide routine formulated the perturbing acceleration which the semidiurnal 

lunar component of the atmospheric tide bulge exerts on the orbit of an artifical satellite. 

The tide bulge was assumed to result from the fact that the earth rotates within the field 

of lunar mass attraction which is inhomogeneous across the terrestrial globe. The routine 

was restricted to the main term of the semidiurnal tide. 

The second tide routine was designed to augment the first by spelling out the 
combined effect on the atmosphere of solar heating and solar mass attraction, plus the 
resulting modification of the earth’s gravitation. The semidiurnal term as well as the diurnal 
term were included. 

The third algorithm was to account for the ocean tide. It started from an existing 
table for the global tide amplitudes and phase angles, and postulated that each tidal height 
value resembles a gravitating point mass that by its presence perturbs the satellite motion. 

Only the semidiurnal lunar tide was considered; however, this algorithm may be regarded as 
setting the pattern for introducing various other spectral components of the fluid tide, if 
desired. 

Further, there was a computer routine for the Newtonian attraction caused by the tide 
bulge of the solid earth. Both the lunar and the solar tide components were included. The 
tidal properties were assumed to be functions of latitude. Accordingly, Love coefficients 
were featured that are zonal harmonic expansions of latitude. An allowance was made for 
tidal lag manifesting itself as a time delay in the response of the tidal mass redistribution 
to the tidal stress field. 
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Several months after completing the coding requests for these four algorithms, an 
opportunity presented itself to devise an interesting alternative to the just mentioned ocean 
tide routine. Namely, a functionalization became available for the tidal ocean surface. This 
consisted of a representation of the tidal ocean surface by an expansion, in surface 
harmonics, that had been least squares fit adjusted to the table of discrete tide amplitudes 
and phases that had been the basis of the ocean tide algorithm. Upon receiving a card 
deck containing the coefficients for this expansion [made available to us through the 
kindness of Dr. C. C. Goad of the National Oceanic and Atmospheric Administration 
(NOAA)], an optional routine was written, which interprets the functionalized ocean surface 
as a continuous, gravitating, tidal mass layer. From it, the perturbing acceleration acting on 
the satellite, due to the presence of the oceanic tide bulge, is then derived via the Poisson 
integral and subsequent gradient formation. While the physics background and the 
mathematical derivations for the original four tide routines are available in the shape of 
formal reports, 5,6,7 ' 8 the background for this latter optional version of the ocean tide 
algorithm had never before been documented. Thus, to preserve it, it is included in this 
report as Appendix A. 


1 - COMPUTER ROUTINE FOR THE PERTURBING ACCELERATION 
CAUSED BY THE LUNAR AIR TIDE 


INPUT DATA 

R = any reasonable value for the mean radius of the earth, in meters 
G = gravitational constant, in meters, kilograms, and seconds 
suggested value: G = 6.6732E - 11 m 3 kg -1 sec -2 
A 2 = a constant associated with the physics of the lunar air tide 
suggested value: A 2 = 0.564 kg nT 2 

The parameters J, n, and t* specify the time instant for which the perturbing acceleration 
due to the tide is to be found. Normally, this is the time instant associated with the 
current time line of the orbit integration. 

J = number of the calendar year 

n = number of the day within the year 

t* = mean solar time at Greenwich (GMT, UTC), in seconds 

x,| inertial, Cartesian components of satellite position vector, in meters, 
x 2 >= associated with the just specified time instant (for a precise definition of 
x 3 J the reference frame, see Appendices B and C of Reference 2) 
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ALGORITHM 


(1 i = j 

60j) = { ^ (1-1) 

lo i 

N = n 6(J, 1975) + (365 + n)6(J, 1976) + (731 + n) 8(J, 1977) + (1096 + n) 6(J, 1978) 

+ (1461+n)5(J, 1979) +(1826+ n)5(J, 1980) (1-2) 

AT = 5.28E - 04 + (3.56E - 08)N (1-3) 

d = 27392.5 +N + AT + t* (M) 

T = d/36525 (1-5) 

s = 270.434358 + 481267.88314IT - 0.001133T 2 + 0.000002T 3- (1-6) 

h = 279.69668 + 36000.768930T + 0.000303T 2 (1-7) 


s and h will result in decimal fractions of degrees. 

Calculate now the earth-fixed, Cartesian satellite coordinates (y p y 2 , and y 3 ) from the 
corresponding inertial coordinates (Xj, x 2 , and x } ). 



= (ABCD)j >n ,t* 



( 1 - 8 ) 


where (ABCD)j n is the transformation matrix that manages the transition from the 

inertial frame of the satellite equations of motion (Basic Inertial System, which is 

associated with either 1950.0 or with the beginning of the day UT of trajectory epoch) to 
the earth-fixed reference frame corresponding to the time instant, J, n, t*. This 
transformation is a computer routine that is already part of TERRA/CELEST. It is external 
to the present algorithm. For details, see Appendices B and C of Reference 2. 
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In the earth-fixed frame, the components of the perturbing acceleration are 


where 


dU 

dr 

3U 

3X 

3U 

30 

3r 

dy, 

' 3X 

dy i 

' 30 

3yi 

a> 

C 

3r 

3U 

3X 

3U 

30 

dr 

3y 2 

' 3X 

9y 2 

' 3 0 

9y 2 

3U 

3r 

3U 

3X 

3U 

30 

3r 

dy 3 

' 3X 

9y 3 

' 30 

dy 3 


(1-9) 


( 1 - 10 ) 


(Ml) 


3r 

yi 

dy i 

r 

3X 

y2 

9yi 

y\ 

3X 

y\ 

dy 2 

+ 

3X 

9y 3 

= 0 

30 

-yi ys 

dy i 

r 2 + y 2 

30 

- yz ys 

dy 2 

r 2 VyJ + y\ 

30 

dy 3 

_ Vyj +_yj 
r 2 


( 1 - 12 ) 

(1-13) 

(M4) 

(1-15) 

(1-16) 

(1-17) 

(1-18) 
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(1-19) 

( 1 - 20 ) 

( 1 - 21 ) 

( 1 - 22 ) 

(1-23) 

(1-24) 

(1-25) 

( 1 - 26 ) 

(1-27) 

(1-28) 

(1-29) 

(1-30) 

(1-31) 

(1-32) 

(1-33) 

(1-34) 

(1-35) 

(1-36) 
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sin X 


_Zl_ 

v'yfTyf 


Finally, transform the perturbing acceleration back into the inertial 



where (ABCD) 1 is the transpose of (ABCD). 


COMPUTER PROGRAM OUTPUT 

The inertial components (T x ] , T Xj , and T Xj ) of the perturbing 
the presence of the lunar air tide are output in terms of m sec' 2 . If 
km sec’ 2 . 


TRIAL DATA FOR PROGRAM CHECKOUT 

R = 6378145 m 

G = 6.6732E - 11 m 3 kg' 1 

A 2 = 0.564 kg m -2 

J = 77 

n = 202 

t* = 50000 sec 

x, = 3151529.23 m 

x 2 = 5458608.75 m 

x 3 = 3639072.50 m 


(1-37) 


frame 


d-38) 


acceleration due to 
required, convert to 


- 2 
sec 
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From external source (CELEST subroutine PRENUT) 

/- 0.8405285753 0.5417623775 0.2289080162E - 02 

(ABCD) = - 0.5417605355 -0.8405316908 0.1413662999E - 02 

\ 0.2689913850E- 02 -0.5190827376E- 04 0.9999963803 

y, = +316648.61 m 

y 2 = - 6290363.38 m 

y 3 = +3647253.32 m 

t** = 208?3333333 

N = (731 + 202) = 933 

AT = 5.612148000E - 04 

t* = 0.5787037037 deg 

d = 2.832607926E + 04 
T = 0.7755257840 
h = 2?819922141E + 04 
s = 3?735060861E + 05 
X = 272?8817616 
h-s = - 3?4530686469E + 05 
a = -3?448331496E + 05 
2a = - 6?896662992E + 05 
x = 0.5011240256 
P^(x) = 2.246624133 
P^(x) = 4.256662025 
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— = 7.070049216E- 12 

dr 


au 

— = - 5.416307419E-04 
a A 


— = 2.46749468IE - 05 
od 


9r 

— = 4.350677408E - 02 

9y, 


— = 1.585715104E - 07 

ayj 

90 

-— = - 3.4615995 18E- 09 

3y i 

3r 

r— = - 8.642811299E- 01 
dy 2 


-— = 7.982281040E- 09 

dy 2 

90 

-— = 6.876619114E- 08 

dy 2 

3r 

r— = 5.011240258E - 01 

9y 3 



90 

r— = 1.189005543E - 07 

ay 3 

T y] = - 8.566502457E - 11 
T y2 = - 8.737156821E- 12 


m sec 


2 - 2 
m sec 


2 - 2 
m sec z 


m 


-1 


m 


-1 


m 


-1 


m 


-1 


m 


-1 


m sec 
ni sec ' 2 
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T y = 6.476836379E - 12 m sec" 2 

T x = 7.675476994E - 11 m sec" 2 

‘ = 7.675476994E - 14 km sec" 2 

T x 2 = - 3.906656638E - 11 m sec" 2 

= - 3.906656638E - 14 km sec" 2 

T XJ = 6.268367431E - 12 m sec” 2 

= 6.268367431E- 15 km sec" 2 

T = |f| = 8.635267078E - 11 m sec" 2 

= 8.635267078E - 14 km sec" 2 


2-COMPUTER ROUTINE FOR THE PERTURBING ACCELERATION 
CAUSED BY THE SOLAR AIR TIDE 


INPUT DATA 


R = any reasonable value for the mean radius of the earth, in meters 
G = gravitational constant, in meters, kilograms, and seconds 
suggested value: G = 6.6732E- 11 m 3 kg" 1 sec" 2 
A, i constants associated with the physics of the solar air tide 
) = suggested values: Aj = 6 kg m" 2 
A 2 j A 2 = 11.9 kg m~ 2 

The parameters J, n, and t* specify the time instant for which the perturbing 
acceleration due to the tide is to be found. Normally, this is the time instant associated 
with the current time line of the orbit integration. 

J = number of the calendar year 
n = number of the day within the year 
t* = mean solar time at Greenwich (GMT, UTC), in seconds 

x, i inertial, Cartesian components of the satellite position vector, in meters, 
x 2 / = associated with the just specified time instant (for a precise definition 
x 3 \ of the reference frame, see Appendices B and C of Reference 2) 
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ALGORITHM 


First, calculate the earth-fixed, Cartesian satellite coordinates (y,, y 2 , and y 3 ) from 
the corresponding inertial coordinates (x r x 2 , and x 3 ). 


= (AElCCDj n , t . ^x 2 ^ (2-1) 

where (ABCD)j n is the transformation matrix that manages the transition from the 

inertial frame of the satellite equations of motion (Basic Inertial System, which is 

associated with either 1950.0 or with the beginning of the day UT of trajectory epoch) to 
the earth-fixed reference frame corresponding to the time instant, J, n, t*. This 
transformation is a computer routine that is already part of TERRA/CELEST. It is external 
to the present algorithm. For details, see Appendices B anc C of Reference 2. 

In the earth-fixed frame, the components of the perturbing acceleration are 



where 


T = ^ 9U 0A_ 8U 80 

y * dr dy, ax 8y, 80 

0U 0r_ + 0U 8A_ 0U 00 

y2 0r 0y 2 0X 0y 2 80 8y 2 

= 0U8r_ + 8U3X_ 8U 80 
J 3 8r 8y 3 8X 8y 3 80 8y 3 


0r_ _ y, 

dy i r 


i = 1,2,3 


8X_ _ y 2 

dy i y] + y\ 


dX = y, 
dy 2 y] + y 2 


( 2 - 2 ) 

(2-3) 

(2-4) 


(2-5) 

( 2 - 6 ) 

(2-7) 
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ax 

0y 3 


— = o 


30 

0y i 

30 

0y 2 

be 

0y 3 

au 

3r 


au 

ax 


au 

30 


- y. y 3 

r 2 V y f + y ^ 

~ y 2 y 3 
r 2\/ y 2 + y 2 

Vy? + y 3 

r 2 

4R4 jr3 

a, -p-Pj(sin0) cos a - a 2 —P 2 (sin 0) cos 2/3 

5r5 

+ a 3 —g- P 2 (sin 0) cos 20 

R 4 2R 3 

a i Pj(sin 0) sin a- a 2 —j- P^(sin 0) sin 20 

2r5 

+ a 3 -jj- P 2 (sin 0) sin 20 

3R 4 , 3R 3 

a, (15 sin^ 0 - 11) sin 0 cos a - a 2 - 3 — sin 20 cos 20 


Pi(x) = 


( 2 - 8 ) 


(2-9) 


( 2 - 10 ) 


( 2 - 11 ) 


( 2 - 12 ) 


(2-13) 


+ a 3 —jj— (7 sin 2 0-4) sin 20 cos 20 

(2-14) 

= + Vyj +y 2 + y] 

(2-15) 

8 jtGR 


= A,- 

1 105 

(2-16) 

5tt 2 GR 


= 

2 64 

(2-17) 

5tt 2 GR a 2 

2 3072 " 48 

(2-18) 

= ~ Vi - x 2 (5x 2 - 1) 

(2-19) 
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P^(x) = 3(1 - X 2 ) 

P 4 (x) = y (1 - x 2 )(7x 2 - 1) 


x 


sin 26 



r 


2 y 3 Vy 2 +y 2 
r 2 


cos a 

sin a 
cos 2 0 
sin 2^ 
cos 0 
sin 0 
cos X 

sin X 


cos a* cos X - sin a* sin X 
cos a* sin X + sin a* cos X 
cos 2 0 - sin 2 0 
2 sin 0 cos 0 

cos 0* cos X - sin 0* sin X 

sin 0* cos X + cos 0* sin X 
Vi 

VyJ + y\ 

y 2 

V y\ + y\ 


a* = t** - 78° 
0* = t**- 146° 


t 


** 


360 ^ 

86400 1 


( 2 - 20 ) 
( 2 - 21 ) 

( 2 - 22 ) 

(2-23) 

(2-24) 

(2-25) 

(2-26) 

(2-27) 

(2-28) 

(2-29) 

(2-30) 

(2-31) 

(2-32) 

(2-33) 

(2-34) 


Finally, transform the perturbing acceleration back into the inertial frame 



(2-35) 
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where (ABCD) T is the transpose of (ABCD). 

While coding the above algorithm for use with TERRA and CELEST, it was found 
helpful to keep in mind the following. A time line is characterized by the time elapsed 
since the epoch (starting instant) of the particular orbital arc. The epoch is specified in 
terms of year, number of the day UT in the year, and time (t EP ), in seconds, elapsed 
from the start of the day UT until epoch time. Usually, t Ep =0. The present algorithm is 
periodic in t*, the latter being the Universal time associated with the individual time line, 
t* starts from zero at the beginning of the day UT during which the time line occurs. It 
is thus unimportant how many integer days have elapsed since epoch. For an illustration, 
see Figure 1. 


EPOCH 


TIME LINE 



EPOCH INTEGER NUMBER OF DAYS DAY UT 

DAY UT 0F 

TIME LINF 


Figure 1. Universal Time in the Algorithm for the 
Perturbing Acceleration Due to the Solar Air Tide 


To actually find the numerical value of t* corresponding to the time line, consider 

that 


t EP = Epoch time (Universal time at epoch) 

At = time difference between time line and epoch 
t* = Universal time of time line 

Now, express t Ep (if non-zero) and At and t* in terms of a common time unit (usually 
mean solar seconds). Find 
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(2-36) 


Express t' in terms of days and decimal fraction of days (1 day = 86400 sec). Now, 

t* = FRACT(t') (2-37) 

Finally, convert t* to angular time 


t** = 360 

which may be substituted for Equation 2-34. 

COMPUTER PROGRAM OUTPUT 

The inertial components (T X] , T X2> and 
the presence of the solar air tide are output 
km sec -2 . 

TRIAL DATA FOR PROGRAM CHECKOUT 

R = 6378145 
G = 6.6732E - 11 
A, = 6 
A 2 = 11.9 
J = 1977 
n = 202 
t* = 50000 


FRACT(t') (2-38) 


T x 3 ) of the perturbing acceleration due to 
in terms of m sec -2 . If required, convert to 


m 

m 3 kg -1 sec -2 
kg m -2 
kg m -2 

sec 
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x, = 3151529.23 


m 


x 2 = 5458608.75 m 

x 3 = 3639072.50 m 


From external source (CELEST subroutine PRENUT) 


(ABCD) = 


0.8405285753 
0.5417605355 
0.2689913850E - 02 


0.5417623775 
“0.8405316908 
"0.5190827376E- 04 


0.2289080162E - 02 
0.1413662999E - 02 
0.9999963803 


*i 

y 2 


y 3 


x 

9 


t** 


0 

Pj(x) 

P\M 

P^(x) 


316648.61 m 

- 6290363.38 m 

3647253.32 m 

7278145.00 m 

272®8817616 

30°07439285 

208°3333333 

403®215095 

335°215095 


0.3318192840 

2.246624133 

4.256662025 


6.112661413E - 04 

m 2 

sec -2 

3.905397704E - 03 

m 2 

sec' 2 

8.136245217E - 05 

m 2 

sec' 2 
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au 

dr 

au 

ax 

au 
a o 


'i 

be 

3y, 

dr 

dy 2 

ax 

ay 2 

3j9_ 

9y 2 

3r 

dy 3 

ax 

dy 3 

be 

dy 3 


= - 1.450823033E- 09 

= 8.799066591E - 03 

= - 6.659222442E - 03 

= 4.350677408E - 02 

= 1.585715 104E - 07 

= - 3.461599518E- 09 

= - 8.642811299E- 01 

= 7.982281040E- 0 9 

= 6.876619116E- 08 

= 5.011240258E- 01 

= 0 

= 1.189005543E - 07 

= 1.355212210E - 09 


t 


m 


-1 


m 


-l 


m 


-l 


m 


-1 


m 


-l 


m sec 


-2 


16 





>2 


T 


8.662262286E - 10 

m sec -2 

- 1.518827519E - 09 

m sec -2 

- 1.612467289E- 09 

m sec -2 

- 1.612467289E - 12 

km sec - 

6.191232115E- 12 

m sec -2 

6.191232115E - 15 

km sec - 

- 1.514495280E - 09 

m sec -2 

- 1.514495280E - 12 

km sec - 

IT 1 = 2.212190101E - 09 

m sec -1 

= 2.212190101E - 12 

km sec - 


3 - FIRST COMPUTER ROUTINE FOR THE PERTURBING ACCELERATION 
CAUSED BY THE M 2 COMPONENT OF THE OCEAN TIDE (BASED ON A 
DISCRETE REPRESENTATION OF THE SEA SURFACE) 


INPUT DATA 

Jv = tidal amplitude on area element ij (these are output data from the 
NSWC/Schwiderski Ocean Tide Program); if available in meters, use 
directly; otherwise, convert to meters 

6 ^ = tidal phase angle on area element ij (these are output data from the 
NSWC/Schwiderski Ocean Tide Program), available in terms of degrees 

NP = number of grid points in the NSWC/Schwiderski Ocean Tide Program; 
this may be expected to be a five-digit integer 

R = radius of the earth (semimajor axis of a suitable reference ellipsoid), 
in kilometers 

suggested value: R = 6378.145 km 
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e 2 = square of eccentricity of reference ellipsoid 
suggested value: e 2 = 0.00669342 

in case it is desired to start from the flattening, f, find e 2 from 
e 2 = (2 - Of 

M e = gravitational constant of the earth, in kilometers and seconds 
suggested value: = 398601 km 3 sec -2 

G = Newton’s constant, in kilometers, kilograms, and seconds 
suggested value: G = 6.6732E-20 km 3 kg " 1 sec " 2 

p = density of water, in kilograms and kilometers 
suggested value: p = l.E+12 kg km " 3 

a = rate of mean mean longitude of moon 

180 . 
suggested value: a = — 1.40519E- 04 deg sec ' 

7T 

The parameters J, n, and t* specify the time instant for which the perturbing 
acceleration due to the tide is to be found. Normally, this is the time instant associated 
with the current time line of the orbit integration 

J = number of the calendar year 
n = number of the day within the year 
t* = mean solar time at Greenwich (GMT, UTC), in seconds 


Xj\ inertial, Cartesian components of satellite position vector, in 
x 2 ( = kilometers, associated with the just specified time instant (for a 
x 3 ( precise definition of the reference frame, see Appendices B and C of 
' Reference 2) 

NMAX = limit on the number of terms in the expansion of the tide potential 


SUBROUTINE FOR THE TIDE POTENTIAL COEFFICIENTS 

This is the first part of the algorithm for the first ocean tide computer routine. The 
time-independent constituents of the expansion coefficients of the ocean tide potential 
(“ F nm > "L’ “ H nm’ and are calculated from the amplitude and phase tables 

produced by the NSWC/Schwiderski Ocean Tide Computer Program. Once established, the 
latter remain unchanged as long as the current versions of the amplitude and phase tables 
remain valid. Thus, the subroutine under discussion will be exercised quite infrequently, 
namely, just once each time a new edition of the just mentioned amplitude and phase 
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tables becomes available. That is expected to occur at intervals of several years, during 
which the present subroutine will remain dormant. 

For the geometry associated with the index ij, see Figures 2, 3, and 4. 


“ij = 1° -3 P G ASjj f y cosSy 


h = 'O^PGAS.j ^ sin 6 y 

- (ii) 2 

= area of the nonpolar surface area element for 
j = 2,3,4 ••• ,j MAX <180 

3 


(AS ij )p OLAR 2 (180) R 


(3-1) 

(3-2) 

(3-3) 


(3-4) 


= area of the polar surface area element (neighboring the North 
Pole) 


(3-5) 

(3-6) 

(3-7) 


0 ii = 

7r 

— (i 

2 

180 \ 2 / 


7T 

/ l\ 

x ij = 

180 

m 

p i , = 

(- 

—sin 2 o\ 

2 ‘V 


To each index ij assign now a counting number v, 1 < v < NP. 


NP 


N P 

T p 2n+1 a P 

H v n nm 

v- 1 

(3-8) 

N P „ 

T p 2n+1 0 P 

Z—i ^ v nm 

v- 1 

(3-9) 


(n “ m)! 1 _ 

-- - V o l 'a h 1 ' 

(n + m)! R 2 " M v v nr 


(3-10) 
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Figure 2. Division of the Earth's Surface Into Area Elements 
According to the Schwiderski Ocean Tide Model 



Figure 3. Position of Point Mass at the Geometrical Center of the 
Surface Area Element Associated With ij 
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V, INDICATES THE POINT WHERE GREENWICH MERIDIAN AND EQUATOR INTERSECT 
y 2 = LOCATED IN THE EQUATORIAL PLANE 


Figure 4. Position of the Point Mass M,j in the 
Earth-Fixed Cartesian Coordinate Frame 




nm 


(n- m)! 1 

(n + m)! M E R 2n 


NP , , 

r p 2n+ '/? h v 

L-». r v r v nm 
v * 1 


(3-11) 


6 


s 

k 


f i e = k 

1 if 

lo t * k 


(3-12) 


Calculate the f*V and h" as follows, noting that 

n m n m ° 


p v = 


_R 

Pv 


(3-13) 


and 


h 1 ' =0 for all values of n 

n o 


(3-14) 
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and that 


= 0 for any n < m 


(3-15) 


Obtain, separately for each v , the required f" m and h^ m from the following 
recurrence relations: 

To advance in n, evaluate 


f "*'[ <2n *' > sin - (" + p, K- .J 6 > 

"n+ i ,m - n - m + 1 [ (2n + » “ 8 . h 'n,m " <" + ■"> P„ H-, ,.] <3-17) 

To advance in m, use 

f n + 1 , n +1 = ( 2n+ Dp^costf, cosX, P „ - cos 0„ sin X, h* J (3-18) 

h n+ 1 ,n+ 1 = (2n+ Dp^cosfl, cos\ v h;_n + cos 0^ sin X, rj (3-19) 

Start these recurrences from 


f" = — 
o.o p 


= R 


sin 0„ 


(3-20) 

(3-21) 


^t ,o 0 


(3-22) 


and terminate the procedure at n = NMAX. 

Store the resulting a F nm , a F nm , “H nm , and a H nm . They will be programmed as 
constant parameters into the subroutine for the perturbing acceleration and will be expected 
to serve for all subsequent runs of that subroutine, until updated. Note that the “F 
^nm' “H nm , anc * are dimensionless quantities. 
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While it would be decidedly impractical to compute P m and h" m from the following 
equation, for program checkout purposes it is possible to find individual, numerical values 
for these quantities from 



R" 


( cos mX y ' 
sin mX^ 


This relationship may also be used to justify Equation 3-15, considering that 

P™(x) = 0 for n < m 

which is one of the basic identities of Legendre function theory. 


(3-23) 


(3-24) 


SUBROUTINE FOR THE PERTURBING ACCELERATION 

First, calculate the earth-fixed, Cartesian satellite coordinates (y,, y 2 , and y 3 ) from 
the corresponding inertial coordinates (x,, x 2 , and x 3 ). 



= (ABCD) 3 n t . 



(3-51)* 


where (ABCD)j n t , is the transformation matrix that manages the transition from the 
inertial frame of the satellite equations of motion (Basic Inertial System, which is 
associated with either 1950.0 or with the beginning of the day UT of trajectory epoch) to 
the earth-fixed reference frame corresponding to the time instant, J, n, t*. This 
transformation is a computer routine that is already part of TERRA/CELEST. It is external 
to the present algorithm. For details, see Appendices B and C of Reference 2. 

Find 


r = + \/x^ + \\ + x^ (3-52a) 

or 

r = +Vyf + y] +y 3 2 (3-52b) 


as convenien* 


t 


Numbering sequence advanced to emphasize parallelism between certain equations of the algorithms in this and in the next chapter. 
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Calculate now the time dependent expansion coefficients for the tide potential from 


F nm = ^ F n m COS K + X) + flF nm sin (at* + X ) (3-53) 

H nm = a H nm cos K + X) +(J H nm sin (at* + X ) (3-54) 


Obtain X as follows 

{ 1 n = m 

if (3-55) 

0 n m 

N v = N 6(J, 1975) + (365 + N)6(J, 1976) + (731 + N) 6(J, 1977) 

+ (1096 + N)6(J, 1978) + (1461 + N) 6(J, 1979) (3-56) 

AT = 5.28E- 04 + (3.56E - 08)1^ 

d 0 = 27392.5 + N l + AT 

X - d o 
0 36525 

X = 270.434358 + 48 1267.88314137 T 0 -.001133 T^+.0000019 T^ (3-60) 

Note that the day number N, and thus X , must be updated whenever t* runs into the 
next day. 

Remember that F nm and H nm are linear functions of two trigonometric terms. To 
evaluate these trigonometric terms, let t* and t* + , be the values of Universal time for 
which subsequent integration steps (time lines) are to be performed. Update the 
trigonometric time factors as follows: 

cos (ot* + , + X ) = cos {(at* + X ) + aAt } = cos (at* + X ) cos aAt - sin (at* + x ) sin oAt (3-61) 


(3-57) 

(3-58) 

(3-59) 
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sin (ot* + , + x) = sin {(at* + x) + oAt } = sin (at* + x) cos oAt + cos (at* + x) sin aAt (3-62) 




Now calculate the Cartesian components of the perturbing acceleration in the 
earth-fixed frame 


N M A X n / at i -\\r 

v V F 3 Au» + H 

*-> Z-f V nm nm a,. 


N M A X n / gu 3Y 

y y (p —A™ + H —*'-$■ 

*—> l—i V nm a v nm a v 

n = 0 m =0 \ ^2 


N M A X n / gjj 


y y F —^ + H — 

Z-i Z-i l nm 3 V nm 

n = 0 m =0 \ °y-\ 


1/1 1 _ 

R \2 ^ fnn ^ n+ l.m-l 2 ^ n+ 1 .m + 1 


1(1 1 

— 1 - — a y -—V 

R \ 2 mn n+1 >m-l 2 n+1 ,m +1 


— (n-m-M)U n+lm 

1 /1 1 N 

- I - A V - - V 

r \^2 m n n+I,m-l j n+l,m + l 

l_ / l + \_ 

r \ 2 A mn n+ l.m-l 2 ^ n+ Cm + l 


— (n- m + 1 )V 

R n + 1 ,r 
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and 


n 

= (n - m + l)(n - m + 2) 

(3-73) 

U n.-. 

(n-D! U nl 

(n + 1)! n 1 n(n + 1) 

(3-74) 

V n,- 1 

(n-1)! V nl 

(n+1)! 1,1 n(n+l) 

(3-75) 


Obtain the values of the individual eigenfunctions from the following recursive 
relationships: 


p = - (3-76) 

r 

sin 0 = — (3-77) 

r 

V =0 for all values of n (3-78) 

n o 

U = V m = 0 for all n < m (3-79) 

n m n m 


To advance in n, evaluate 

U . . = --- j(2n + 1) sin U - (n + m)p U . 1 

n+l,m n-m+1 < 7 nm ' n-l,mi 

V . = --- | (2n + I) sin \j/ V - (n + m)p V . \ 

n+l,m n-m+1 ' nm ' n-l,m> 

To advance in m, use 

- <2n+ Dp(7 U„- (3-82) 


(3-80) 

(3-81) 


26 





(3-83) 


- «n + OP^ V„„ + ^ U„„) 

Start from 


Vo 

Me 

r 

(3-84) 

U .,o 

R y 3 

M e f3 

(3-85) 

Vo 

ii 

< 

o 

II 

O 

(3-86) 


The Cartesian components of the perturbing acceleration, in the earth-fixed reference 
frame, are now 


• _ 

3 0 

(3-87) 

y i 

3yj 

y 2 

30 

ay 2 

(3-88) 

> _ 

30 

(3-89) 

y 3 

9y 3 


Finally, rotate the perturbing acceleration back into the inertial coordinate system 


= ( ABCD >J,n,t. 



(3-90) 


where (ABCD) T is the transpose of (ABCD). 
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COMPUTER PROGRAM OUTPUT 


The inertial components (T x ^, T Xj , and T x ^) of the perturbing acceleration due to 
the presence of the ocean tide are output in terms of km sec' 2 . If required, convert to 
m sec -2 . 


TRIAL DATA FOR PROGRAM CHECKOUT 

Non-zero tidal heights and associated phase values were assumed for nine surface area • 
elements. Several of the “F nm , ^F , °H nm , and w H nm were subsequently computed. The 
algorithms for the perturbing acceleration, featured in the first and second computer 
routines for the ocean tide, are identical. Trial data for these segments of the two tide 
routines will therefore be listed only once, namely, for the second ocean tide routine, 
below. 






Quantities Associated With the Surface Area Elements 




r 

t> 

e 

A 

P 


a 

V 

K 

><i 

» 

m 

deg 

rad 

rad 

km 

km J 

km J sec 1 

km’ sec' 1 

u 

1 

10 

25 

1.562069681 

8.7266463 E 03 

6356.800824 

108.1411251 

6.5403462E - 08 

3.0498135E - 08 

2,1 

2 

10 

25 

1.562069681 

2.61799388L 02 

6356.800824 

108.14 1 1251 

6.5403462E - 08 

3.0498135E - 08 

3.1 

3 

10 

25 

1.562069681 

4.36332313Ii - 02 

6356.800824 

108.14 11251 

6.5403462E - 08 

3.0498135E - 08 

1.2 

4 

to 

25 

1.544616388 

8.7266463 E. - 03 

6356.813825 

432 4766612 

2.6156072E - 07 

1.2196 7 7 7E - 07 

2,2 

5 

20 

30 

1 544616388 

2.61799388E - 02 

6356.813825 

432.4766612 

4.9987043E - 07 

2.8860033E - 07 

3.2 

6 

20 

30 

1.544616388 

4.36332312b - 02 

6356.813825 

432.4766612 

4.9987043E - 07 

2.8860033E - 07 

1.3 

7 

10 

25 

1.527163095 

8.7266463 E - 03 

6356.839812 

648.550316 

3.9224149E - 07 

1.8290521E - 07 

2.3 

8 

20 

30 

1.527163095 

2.61 799388E 02 

6356.839812 

648 550316 

7.496153 E - 07 

4.327906 E - 07 

3,3 

9 

20 

30 

1.527163095 

4.36332313E - 02 

6356.839812 

648.550316 

7.496153 E - 07 

4.327906 E - 07 
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Values for the 

f and h^ m 



f' 

‘o,o 

h o.o 

f l,0 

h i ,0 

V 

km 1 

km 1 

km" 1 

km 1 


1 

1.5731184E-04 

0 

1.5783403E - 04 

0 

2 

1.5731 184E - 04 

0 

1.5783403E - 04 

0 

3 

1.5731 184E -04 

0 

1.5783403E - 04 

0 

4 

1.5731 15 IE - 04 

0 

1.5778531E- 04 

0 

5 

1.5731 15 IE-04 

0 

1.577853IE - 04 

0 

6 

1.5731 15IE-04 

0 

1.5778531E — 04 

0 

7 

1.5731087E - 04 

0 

1.5768788E-04 

0 

8 

1.5731087E - 04 

0 

I.5768788E -04 

0 

9 

1.5731087E - 04 

0 

1.5768788E - 04 

0 


V 

P 

i.i 

km' 1 


h u 
km 1 


f 2,0 
km 1 


h 2 ,o 

km 1 

1 

1 3773442E 

-06 

1.201990 IE 

-08 

1.5835193E 

-04 

0 

2 

1.3769246E 

-06 

3.605604 IE 

-08 

1.5835193E 

-04 

0 

3 

1.3760857E 

-06 

6.0081198E 

-08 

1.5835193E 

-04 

0 

4 

4.1315963E 

-06 

3.6055895E 

-08 

1.5820627E 

-04 

0 

5 

4.1303378E 

-06 

1.0815670E 

-07 

I.5820627E 

-04 

0 

6 

4.127821 IE 

-06 

1.802 24 56E 

-07 

1.5820627E 

-04 

0 

7 

6.8845393E 

-06 

6.0080464E 

-08 

I.57915I3E 

-04 

0 

8 

6.8824422E 

-06 

1.8022309E 

-07 

1.5791513E 

-04 

0 

9 

6.8782468E 

-06 

3.0031082E 

-07 

1.5791513E 

-04 

0 



f 2,1 
km' 1 


h 2,l 

km 1 


f 2.2 
km 1 


h 2.2 
km 1 


1 

4.1457488E 

-06 

3.6179402E 

-08 

3.61752671: 

-08 

6.3 144 1631: 

- 10 

2 

4.1444859E 

-06 

1.0852718E 

-07 

3.6131 193E 

-08 

1.8935 556E 

-09 

3 

4.1419607E 

-06 

1.8084191E 

-07 

3.6043099E 

-08 

3.1533625E 

-09 

4 

1.24321 20E 

-05 

1.0849347E 

-07 

3.2550933E 

-07 

5.68 1 78641: 

-09 

5 

I.2428333E 

-05 

3.2544735E 

-07 

3.251 1274E 

-07 

1.70384371: 

-08 

6 

1.2420760E 

-05 

5.4230210E 

-07 

3.24320061: 

-07 

2.8374329E 

-08 

7 

2.07031I6E 

-05 

1.80673 35 E 

-07 

9.038143 IE 

-07 

1.57761371: 

-08 

8 

2.0696809E 

-05 

5.4I96503E 

-07 

9.0271315E 

-07 

4.730919 IE 

-08 

9 

2.0684198E 

-05 

9.030916 IE 

-07 

9.0051217E 

-07 

7.8 784 606E 

-o.s 


Values for the f 


and h" (Continued) 


V 

fK 

'3.0 

km' 1 

h 3,0 
km - 1 

f 3,l 
km -1 

^3,1 
km - 1 

1 

1.5886548E-04 

0 

8.3188628E - 06 

7.2597615E-08 


I.5886548E-04 

0 

8.3163287E - 06 

2.1777073E - 07 

3 

1.5886548E-04 

0 

G.3112614E-06 

3.628775 IE-07 

4 

1.5857391E - 04 

0 

2.4934851 E-05 

2.1760315E - 07 

5 

1.5857391E - 04 

0 

2.4927255E - 05 

6.5274316E - 07 

6 

1.5857391E - 04 

0 

2.4912067E - 05 

1.0876843E - 06 

7 

1.5799154E - 04 

0 

4.1485684E - 05 

3.6204008E - 07 

8 

1.5799154E - 04 

0 

4.1473047E - 05 

1.0860100E -06 

9 

1.5799154E - 04 

0 

4.1447777E - 05 

1.8096490E -06 


V 

f 3 ,2 
knr 1 

h 3,2 

knT 1 

ev 

3,3 

km' 1 

h 3.3 
km - 1 

1 

1.8147675E - 07 

3.1676885E - 09 

1.5834220E - 09 

4.1463365E- 11 

2 

1.8125565E - 07 

9.4992060E - 09 

1.579082 E - 09 

1.2427645E- 10 

3 

1.8081372E - 07 

1.5819151E - 08 

1.5704138E - 09 

2.0674889E - 10 

4 

1.6324485E - 06 

2.8494495E -08 

4.2739029E - 08 

1.1191609E - 09 

5 

1.6304596E - 06 

8.5448767E - 08 

4.2621885E -08 

3.3544151E - 09 

6 

1.6264843E - 06 

1.4229893E - 07 

4.2387916E- 08 

5.58047 E - 09 

7 

4.5299018E - 06 

7.9069730E - 08 

1 9774213E- 07 

5.1780599E-09 

8 

4.5243828E - 06 

2.371 1 286E — 07 

1.97200I3E - 07 

1.55 19987E - 08 

9 

4.5133516E - 06 

3.9486710E-07 

1.9611762E-07 

2.5819375E - 08 


V 

'4,3 

km' * 

h 4 ,3 
km' 1 

1 

l.l I20747E-08 

2.9120701E - 10 

2 

I.I090266E-08 

8.7282284E - 10 

3 

1.1029387E - 08 

1.4520463E -09 

4 

3.0007426E - 07 

7.857721 IE-09 

5 

2.9925178E -07 

2.3551626E - 08 

6 

2.9760907E - 07 

3.9180977E - 08 

7 

1.3875 122E -06 

3.6333286E - 08 

8 

1.383709IE - 06 

1.0890027E - 07 

9 

1.3761134E-06 

1.8116877E -07 
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Time-Independent Coefficients of the Tide Potential 


n 

m 

“F 

n m 


"F 

n m 


“H 

n m 


"H 

n m 


0 

0 

8.4018456F. - 

12 

4.6140106E — 

12 

0 


0 


1 

0 

8.368164 IE — 

12 

4.5954868E - 

12 

0 


0 


I 

1 

2.9297877F, — 

13 

I.6202500E- 

13 

4.3123569E - 

15 

2.4544968E - 

15 

2 

0 

8.3290386E - 

12 

4.5739464E — 

12 

0 


0 


2 

I 

2.9177584E - 

13 

1.6I35948E — 

13 

4.2946458E - 

15 

2.4444118E- 

15 

2 

2 

2.78405 15E- 

15 

1.5423222E- 

15 

8.2132887E - 

17 

4.683368 E- 

17 

3 

0 

8.2845357E - 

12 

4.5494 263E- 

12 

0 


0 


3 

1 

2.9046632E - 

13 

1.6063487E - 

13 

4.2753632E - 

15 

2.4334303E - 

15 

3 

2 

2.7724443E - 

15 

1.5358913E - 

15 

8.1790436E — 

17 

4.6638388E — 

17 

3 

3 

1.5813081E - 

17 

1.0259185E - 

17 

8. 206853 5E - 

19 

4.6816230E - 

19 

4 

3 

1.8435 105E- 

17 

1.0215973E- 

17 

8.1722863E- 

19 

4.6619036E - 

19 


4-SECOND COMPUTER ROUTINE FOR THE PERTURBING ACCELERATION 
CAUSED BY THE M 2 COMPONENT OF THE OCEAN TIDE (BASED ON A 
FUNCTIONALIZATION OF THE SEA SURFACE) 


INPUT DATA 


NOAA/Goad expansion coefficients of the tidal surface, in meters; 
these are specified by a computer card deck 

R = radius of the earth (semimajor axis of a suitable reference ellipsoid) 
in kilometers 

suggested value: R = 6378.145 km 

G = Newton’s constant, in kilometers, kilograms, and seconds 
suggested value: G = 6.6732F - 20 km 3 kg 1 sec 2 

p ( = gravitational constant of the earth, in kilometers and seconds 
suggested value: p t = 398601 km 3 sec 2 

p = density of water, in kilograms and kilometers 
suggested value: p = l.E+12 kg km 3 


Ci, 

S ., 

S u 
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o = rate of mean mean longitude of moon 
180 

suggested value: u = — 1.40519k 04 deg see 1 

7T 

NMAX = limit on the number of terms in the expansion of the tide potential; 

NMAX is equal to the largest value of the index i for the expansion 
coefficients of the tidal surface 

The parameters J, n, and t* specify the time instant for which the perturbing 
acceleration due to the tide is to be found. Normally, this is the time instant associated 
with the current time line of the orbit integration. 

J = number of the calendar year 
n = number of the day within the year 
t* = mean solar time at Greenwich (GMT, UTC), in seconds 

X|| inertial, Cartesian components of satellite position vector, in 

x 2 l = kilometers, associated with the just specified time instant (for a 
x 3 I precise definition of the reference frame, see Appendices B and C of 
' Reference 2) 


SUBROUTINE TOR THE TIDE POTENTIAL COEFFICIENTS 


This is the first part of the algorithm for the second ocean tide computer routine. 
The time-independent constituents of the expansion coefficients of the ocean tide potential 
(F^, F”, Hjj, II'') are calculated from the expansion coefficients for the NOAA/Goad tidal 
ocean surface. Once established, the latter remain unchanged as long as the current version 
ol the NSWC/Schwiderski Ocean Tide Model (on which the NOAA/Goad tidal surface is 
based) remains valid. Thus, the subroutine under discussion will be exercised quite 
infrequently, namely, just once each time a new edition of the NSWC/Schwiderski Ocean 
Tide Model becomes available. That is expected to occur at intervals of several years, 
during which the present subroutine will remain dormant. 


u 


u 


2ttR j Gp 

2 

Me 

2i + 1 

27rK“Gp 

2 

Me 

2i + 1 

27tR-g fj 

■> 

Ml 

2i + 1 


K) 

10 jq, 
10 \s 

IJ 


(4-1) 

(4-2) 


II' 


(4-3) 



1 


h:: 


27rR 2 Gp 


2i + 1 


10 3 S 


(4-4) 


SUBROUTINE FOR THE PERTURBING ACCELERATION 

Proceed exactly as in Equations 3-51 through 3-90 from the first ocean tide computer 
routine, but delete Equations 3-53 and 3-54 and replace them with 


= F nm cos < at * + X) + F'' m sin (at* + X) 


(4-53) 


H = H' cos (at* + x) + H” sin (at* + x) (4-54) 

n in nm v i nm ' i ' ' 

COMPUTER PROGRAM OUTPUT 

The inertial components (T x ( , T X2 , and T Xj ) of the perturbing acceleration due to 
the presence of the ocean tide are output in terms of km sec -2 . If required, convert to 
m sec -2 . 


TRIAL DATA FOR PROGRAM CHECKOUT 

Non-zero values were assigned to selected NOAA/Goad expansion coefficients, and the 
corresponding expansion coefficients of the tide potential were computed. The algorithm for 
the perturbing acceleration was then exercised with the following results. 


J 

= 1977 



n 

= 202 



t* 

= 50000 

sec 


R 

= 6378.145 

km 


G 

= 6.6732E - 20 

km 3 

kg -1 

H 

= 398601 

km 3 

sec - 2 
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p = 10 


kg km 3 

a = -1.40519E-04 deg sec"" 1 

7r 

NMAX = 4 

C 20 = + 0.2906060089E-01 m 
C 40 = - 0.107121752E + 00 m 
C 43 = + 0.435761219E - 04 m 

C 2 ' 0 = - 0.4424413130E- 01 m 
C 4 ' 0 = + 0.873468034E - 01 m 
C 43 = - 0.160563906E-02 m 

s 20 - o 

S 40 ~ 0 

S 43 = - 0.363303008E - 02 m 

S 20 = 0 
S 40 = 0 

S 4 ' 3 = - 0.264356490E - 02 m 

F 20 = + 4.9742658E - 10 
F 40 = - 1.0186607E - 09 
F^ 3 = + 4.1438160E - 13 
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F 20 = ~ 7.573211 HE - 10 
F" 0 = + 8.30613340E- 10 
F" 3 = - 1.526862IE - 11 

H' 0 = 0 

h; 0 = 0 

= - 3.4547839E - 11 

h 20 = o 

H 40 = 0 

H” 3 = - 2.5138645E - 11 

N t . = 933 
AT = 0.000561 
d Q = 28325.50056 
T 0 = 0.7755099401 
X = 373498.4609 
a = 0.0080511456 
ot* = 402.557282 
cos(ot* + x) = -0.754501 
sin (at* + x) = - 0.656298 

F 20 = + 1.2171968E- 10 

h 20 = 0 

F 40 = + 2.2345060E - 10 


deg/sec 

deg 


35 






0 


From external 


(ABCD) = 


F 43 = + 9.7081216E - 12 
H 43 = + 4.2564846F - 11 

km 
km 
km 


x, = 3151.52923 
x 2 = 5458.60875 
x 3 = 3639.07250 


-0.8405285753 
-0.5417605355 
0.2689913850E- 02 


0.5417623775 

- 0.8405316908 

- 0.5190827376E - 04 


0.2289080162E - 02 
0.1413662999E- 02 
0.9999963803 


y, = + 316.64861 km 

y 2 = - 6290.36338 km 

y 3 = + 3647.25332 km 


9U 2U 

r —— = - 0.0000964047 km/sec 2 

ay, 

av 2n 

r- - - 4.7E - 13 km/sec 2 

ay, 

au 20 

t —~ = +0.0019151219 km/sec 2 

ay 2 

av 2n 

r- = - 2.4E - 14 km/sec 2 

ay 2 
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\0 \Q 


Values for the Eigenfunctions of the Tide 
Potential, in km 2 sec' 2 


n 


m 


U 


n m 


v 


n m 


0 

-1 


0 

0 

54 

1 

1 

-1 

0 

-1 

24 

1 

1 

2 

2 

-1 

-0 

2 

0 

-5 

2 

1 

2 

2 

2 

-94 

3 

-1 

-0 

3 

0 

-16 

3 

1 

0 

3 

2 

-206 

3 

3 

-53 

4 

-1 

0 

4 

0 

-9 

4 

1 

-2 

4 

2 

-136 

4 

3 

-165 

4 

4 

1863 


o o 


.76683963 

0 

.044042677 

-20.74036524 

.05119113 

0 

.088085354 

-41.48073047 

.4584976993 

-9.108257692 

.186455141 

0 

.750986196 

-54.64954615 

.01442066 

-9.489169689 

.0512402699 

-1.017910413 

.10992273 

0 

.6148832385 

-12.21492495 

.435027 

-20.83613333 

.85812191 

354.226451 

.109342534 

2.172137348 

.393545781 

0 

.18685068 

43.44274695 

.7982658 

-13.80747707 

.56486 

1088.924921 

.67849 

380.085929 


5 

5 

5 

5 

5 

5 


-1 

0 

1 

2 

3 

4 


0.0917032956 

2.472201758 

-2.751098868 

136.846714 

-182.4236505 

7366.011828 


1.821726148 

0 

54.65178443 
13.8123671 
1 199.805712 
1502.253452 


-1 

0 

1 

2 

3 

4 


0.0153004451 
8.002095405 
-0.6426186943 
349.1179421 
45.32033435 
11350.89177 


0.303950046 

0 

12.76590193 

35.23756643 

-298.0731712 

2314.94556 
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9Uh, 

9y 3 

aV 20 

9y 3 

9^40 

9y, 


9V 4 q 
9y i 

9U 4 q 

9y 2 

9V 4 q 

9y 2 


9U 4 q 

9y 3 

9 V 4 o 

9y 3 

9U 43 
9y i 


9V 4 3 
9y i 

9U 43 

9y, 


= +0.0075774019 

= 0 

= +0.0004313321 

= + 7.8E- 13 

= - 0.0085686018 

= 0 

= - 0.0019380257 

= 0 

= - 0.5130748473 

= - 0.1112689700 

= - 0.1242624348 


9V 43 

= + 0.6418082461 
9y 2 


km/sec 2 


km/sec 2 


km/sec 2 


km/sec 2 


km/sec 2 


km/sec 2 


km/sec 2 


km/sec 2 


km/sec 2 
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3U 4 3 

-— = +0.0572027292 km/sec 2 

ay 3 

3 V 4 3 

*- = - 0.3762240313 km/sec 2 

dy 3 

T y = - 9.632495E- 12 km/sec 2 

T y2 = + 2.443056E - 11 km/s^c 2 

T yj = - 1.496932 IE - 11 km/sec 2 

T x = - 5.179392E - 12 km/sec 2 

T X2 = - 2.5752406E- 11 km/sec 2 

T X3 = -1.495676E- 11 km/sec 2 


5 - COMPUTER ROUTINE FOR THE PERTURBING ACCELERATION 
CAUSED BY THE TIDES OF THE SOLID EARTH 


INPUT DATA 

The parameters J, n, and t* specify the time instant t for which the perturbing 
acceleration due to the tide is to be found. Normally, this is the time instant associated 
with the current time line of the orbit integration. 

J = number of the calendar year 
n = number of the day within the year 
t* = mean solar time at Greenwich (GMT, UTC), in seconds 

At = tidal lag parameter, in seconds 

w = earth’s sidereal rate of rotation 

suggested value: to = 4.178074622E - 03 deg sec -1 
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the inertial, Cartesian coordinates of the moon, in kilometers, for the 
time instant (t - At), from the TERRA/CELEST moon ephemeris; the 
reference frame is relerred to mean equator and mean equinox at time 
instant t, 

the inertial, Cartesian coordinates of the sun, in kilometers, for the time 
instant (t - At), trom the TERRA/CELEST sun ephemeris; the reference 
frame is relerred to mean equator and mean equinox at time instant tj 

epoch time for the moon ephemeris and sun ephemeris, in Julian Date 
t, = JD2433282.423 

inertial. Cartesian components of satellite position vector, in kilometers, 
associated with the time instant t(J, n, t*) (for a precise definition of 
the reference frame, see Appendices B and C of Reference 2) 


Cartesian satellite velocity components associated with time instant t, in 
km sec' 1 


gravitational constant ot the earth, in kilometers and seconds 
suggested value: /a F = 398601 km 3 sec' 2 

radius of the earth (semimajor axis of a suitable reference ellipsoid) in 
kilometers 

suggested value: R = 6378.145 km 

square of eccentricity of reference ellipsoid 
suggested value: e 2 = .006693421623 

in case it is desired to start from the tlattening, f, find e 2 from 
e 2 = (2 - Of 

expansion coefficients for the Love numbers 
suggested values: k 2() = 0.3 
k 2 , = 0.01 
k 30 - 0.1 
k 3| = 0.01 
k 22 = 0.1 
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m' = mass of moon or sun, in kilograms 

suggested values: m ' L = mass of moon = 7.369328IE + 22 kg 
m^ = mass of sun = 1.99E + 30 kg 

M = mass of earth, in kilograms 

suggested value: M = 5.9731613E + 24 kg 

a' = semimajor axis of moon or sun orbit relative to earth, in kilometers; 
select values that are compatible with the moon ephemeris or sun 
ephemeris (optionally, values may be obtained from a current edition of 
the AMERICAN EPHEMERIS ) 
suggested values: a ' L = 384400 km 

a' s = 149500000 km 


SUBROUTINE FOR THE POSITION OF THE FICTITIOUS MOON 





(5-1) 


cos £ 0 cos 0 cos z - sin £ 0 sin z 

(5-2) 

- sin £ 0 cos 6 cos z - cos £ 0 sin z 

(5-3) 

- sin 0 cos z 

(5-4) 

cos £ 0 cos 0 sin z + sin £ 0 cos z 

(5-5) 

- sin £ () cos 0 sin z + cos £ 0 cos z 

(5-6) 

- sin 0 sin z 

(5-7) 

cos £ 0 sin 6 

(5-8) 

- sin £ 0 sin 0 

(5-9) 

cos 6 

(5-10) 
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(5-11) 


( n = (2304.250 + 1.396T n )T + 0.302T 2 + 0.018T 3 - 

0 ' 0 ’ 3600 


2 ' *» + !°- 79it! S5^ 


6 = j (2004.682 - 0.853T 0 )T- 0.426T 2 - 0.042T 3 j 


3600 


(5-12) 

(5-13) 


£ 0 , z, and 0 result in degrees. 


(t, " t, ) 


= (t, - 2415020.313)(0.273790926497E - 04) 

(5-14) 

= (0.273790926497E - 04)(t, - t,) 

(5-15) 

= ) (J - 50)365 + A 0 + (n - 1) + .077 | 

(5-16) 

IX ,/J-53 \ 


= ,NT br + j 

(5-17) 


The Cartesian coordinates of the fictitious moon (x*, y*, and z*) referred to the 
coordinate system in which the satellite orbit integration is performed may rtow be 
computed: 


- 


( x*\ /coscoAt 
y* I = ( s ‘ n <*>At 
z*/ \ 0 


-sin coAt 
cos o)At 
0 



(5-18) 


Also, compute 


r* = + v/(x*) 2 + (y* ) 2 + (z* ) 2 



z 


* 

I. 


r 


* 

L 


(5-19) 

(5-20) 

(5-21) 

(5-22) 
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SUBROUTINE FOR THE POSITION OF THE FICTITIOUS SUN 





(5-23) 


Evaluate the matrix elements X x , Y x , etc., according to Equations 5-2 through 5-17. 
Then, calculate the Cartesian coordinates of the fictitious sun (x£, y*, z*) referred to the 
coordinate system in which the satellite orbit integration is performed: 



cos toAt 
sin cuAt 
0 


-sin coAt 0\ /x's \ 

coscuAt 0 )( y' s ) 

0 1/\z' s / 


(5-24) 


Further, compute 


+ >/(x*)2 +(y*)2 + (z*)2 



SUBROUTINE FOR SATELLITE POSITION 

r = +x/xf + x^ + x^ 



(5-25) 

(5-26) 

(5-27) 

(5-28) 


(5-29) 

(5-30) 
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■2 


(5-31) 


li 2 — 1 0 )2 4. \2 . \2 


(:■£) 


SUBROUTINE FOR THE PERTURBING ACCELERATION 

T x J , T X2 , and T Xj are the Cartesian components of the perturbing acceleration, in 
km sec- 2 . These are referred to the (inertial) reference frame in which the satellite orbit 
integration is performed, and are thus readily suited for insertion into the TERRA/CELEST 
equations of motion. 


T x, = 71 [ c ° + + TT (5V 2- 2P ^ + 7( 7V 3- F)+ 7< 9V 4-H)] x , 

[c, , c 2 c 3 c 4 1 I 

- [y A 3 + 7^11 + 7^12 + r yPi 3 J j (5-35 

-Mf 3C, V, c, C 3 c, 1 

T x 2 = 7 j [ C o + — + r T(5V 2 - 2Pq)+ r -f(7V 3 -F) + ^(9V 4 -H) x 2 

c, , c 2 c 3 c 4 "I / 

■ 7 A 4 + 7 P 2 1 + 7^22 + 7^23 J j (5-36 

T .J - ^ ■ II, ,J 


-ff 


P 3 1 + r 2 ^3 2 + r 3 p 3 
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Pyx 

= 2XP' + /liP' + yPj 

(5-38) 

Pxj 

= (1 - 5i> 2 )S' + 3(X 2 - n 2 )S' A + 2X(3a«s; + i>S' 6 ) + ni>S' 7 

(5-39) 

Px J 

= v (1 - 7 - 1 ' 2 ) r; + 3wx 2 - n 2 >t; + 6x*ii/r; + (1 - iv 1 >( 2 xt; + mt; » 

(5-40) 

P 21 

= -2/iP' + XP' 2 + vP' 

(5-41) 

P 22 

= a - 5v 2 >s; - 6x^s; + 3<x 2 - m 2 is; + ^xs; - 2 ms;, > 

(5-42) 

Pa 

= */» - n -t^j r; - bXuvV* + 3wx 2 - m 2 it; + (I - iv l mxt; - 2 mt; > 

(5-43) 

Pm 

= -6< +xp; +< 

(5-44) 

Pm 

= lOi'lXS’, + n$' 2 ) + (3 - 1 5 m 2 (S', + (X 2 - fj 2 )S^ + XmS^ 

(5-45) 

Pm 

= (1 - lv l M XT' + juT' 2 ) - 20m( 3 - Iv 1 IT', 

+ X(X 2 - 3m 2 (T^ + /j(3X 2 - m 2 )T',, 

- i4i^((x 2 - m 2 >t; +x/jt;i 

(5-46) 

V 

= 2 <xs; +ms; + 3ms;> 

(5-47) 

II 

= 2 [i'(XT' + //Tj ( + 6< 1 - Si' 2 it; + (X 2 - /j 2 )T^ + X/iT^] 

(5-48) 

where 



<0 

^3 e ^2 0 ^2 2^ uKA o 

(5-49) 


= Kak ; , a 2 

(5-50) 


= Ka 2 a J 

(5-51) 


= Ka J a 4 

(5-52) 


= Ka 4 a 5 

(5-53) 
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and 


k 


— n 2 a 2 a' 3 d * 3 
M 


a J n 2 = p = M(i 
, K 

Ot = — 

a 

K 

a = — 


and 


<** = 


and 


V i = ^ (A ',* +A > 4A >> 

V = 1’ P' P 

2 il o " " 

V = I S' s 

n 1 n n 

V = V x'T 

4 Ml" " 

A' = ~(I if * 1 ) 

0 4 

A' = - (X* 2 /j* 2 ) 

1 4 

A' 2 = 3X*/j* 

A' ( = 3X*i>* 

A^ = in*v* 


(5-54) 

(5-55) 

(5-56) 

(5-57) 

(5-58) 

(5-59) 

(5-60) 

(5-61) 

(5-62) 

(5-63) 

(5-64) 

(5-65) 

(5-66) 

(5-67) 

(5-68) 
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and 


B', = ~ A*( 1 - 5v* 2 ) (5-69) 

n 

, 3 

B 2 = “#i*M - 5v * 2 ) (5-70) 

O 

B 3 = 7 "*(3 - 5^* 2 ) (5-71) 

B 4 = ^A*(A* 2 -3 m* 2 ) (5-72) 

O 

B 5 = “M*(3X* 2 - tx* 2 ) (5-73) 

O 

B^ = ~ v*(\* 2 - n * 2 ) (5-74) 

B, = 15(5-75) 

and 

r ; • ( k ., + ? k . 1 )( i - i ei ) A i + ,s - 76 ' 

P 0 = 1 - 3p 2 (5-77) 

P '. = ( k 2o- ^ k 2 2 ) (’ " ^ 2 )a', + ^ 3 , a '^ B l (5-78) 

P, = A 2 - #J 2 (5-79) 

p ; = ( k 2 0 - ^ k 2 2) (^ ~ a; + ^k 31 a'(?*B; ( 5 - 80 ) 

P 2 = A/i (5-81) 

p ; = ( k 2o + 7 k 2 2 )( | - H^a;- (5-82) 

Pj * Av (5-83) 

k - ( k ,» * f k >0(' - £•')*; • “ k ,,“'u‘ B ; <5-84) 

P 4 = (5-85) 




47 




aiul 


S > = " " k 2. A J +k 30 0, '^ B ', (5-86) 

S, = Ml - 5 v 2 ) (5-87) 

S 2 = - 5 - k 2l A 4 +k 30^* B 2 (5-88) 

5 2 = n( 1 - 5 v 2 ) (5-89) 

5 3 = | k 2 I A 0 +k 3 0 CV '( 3 * B 3 (5-90) 

5 3 = v(3 - 5v 2 ) (5-91) 

s; = k 3O a'0*B; ( 5 . 92 ) 

5 4 = X(X 2 - 3m 2 ) (5-93) 

s ; = k 30 a'(J*B; (5-94) 

S s = /a(3X 2 - /a 2 ) (5-95) 

S 6 = k 2i A '. +k 3o“'( 3 * B l (5-96) 

5 6 = i>(X 2 - n 2 ) (5-97) 

5 7 = k 2 i A 2 + k 3o a '0* B 7 (5-98) 

S 7 = X^ (5-99) 
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and 


t; 

T , 

t; 

t 2 

t; 

T 3 

n 

T 4 

t ; 

T s 

Tl 


1 15 ,/ 1 \ 9 1 , 15 

H e ( k - + 7 k -r M k -J A 3 + T k - a ^ B - 


\i> 1 - - v 2 

3 


15 


( k 20 + 7 k 2’j ]4 k 2 2^ A 4 + 1 k 3| a0*B 2 

^(i- r 2 ) 

| 1 4 ^2 0 + 7 k 2 2^ e 7Q k 2 :| A 0 ' 7 k 3 1 01 

3 - 30v 2 + 35t> 4 

k 31 «'/3*B; 

\v(\ 2 - 3/i 2 ) 

k 3, a '^ B S 

mK3X 2 - m 2 ) 

14 ( k 2 0 ‘ 7 k 22) e ~ 4 k 2 2 1 A 1 ~ 7 k 3i a ^* B 6 


r 6 = tx 2 - n 2 )(i - iv 2 ) 

1' 


!-e 2 ! 


-k ) 

| 14 ' 

( k2 0 ~ 7 kj2 j 

- 4 k 22j 


T, = Xp( 1 - lv 2 ) 


(5-100) 

(5-101) 

(50102) 

(5-103) 

(5-104) 

(5-105) 

(5-106) 

(5-107) 

(5-108) 

(5-109) 

(5-110) 

(5-111) 

(5-112) 

(5-113) 


Note that the components T v of the perturbing acceleration depend, essentially, on 
the parameters 


X*, n*, v*, r*. a', m', x,, x 2 , x 3 , 
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T x . = T x .(X*, m*. v*, r*. a', m', x,, x 2 , x 3 ) 


To obtain the perturbing acceleration due to the moon tide, substitute for these parameters 
those specifying lunar position, orbital semimajor axis, and mass. 


(T x .) l = T x .(X*, n *, v*, r*, m' L , x,, x 2 , x 3 ) 


Otherwise, substitute the solar parameters to produce the perturbing acceleration caused by 
the tidal action of the sun. 


(T x .) s = T x .(\*, p*, n*, r|, a^, m^, x,, x 2 , x 3 ) 

COMPUTER PROGRAM OUTPUT 

The inertial components (T X( ) L , (T x 7 ) L , (T x 3 ) L , (T X] ) S , (T Xi ) $ , and (T X3 > S of the 
perturbing acceleration due to lunar and solar tides of the solid earth are input in terms of 
km sec -2 . 


TRIAL DATA FOR PROGRAM CHECKOUT 

The following trial case is restricted to the lunar tide. While executing this case, the 
entire algorithm was exercised and it was, thus, unnecessary to also consider the solar tide. 


Input Data 

t ••• J = 1977 
n = 88 

t* = 57000 sec 

At = 100 sec 

£3 = 4.178074622E - 03 deg sec -1 
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Epoch for the moon ephenieris (reference frame 2j) 
t, = 1950.0 = JD2433282.423 


(1950.0 means beginning of Besselian year 1950) 

Epoch for the satellite orbit (reference 

frame 2 ? ) 

V 

•• J = 1977 



n = 88 



t* = 0 

sec 

*2 

= JD2443231.5 


Moon coordinates, in frame 2., at time 

: (t- At) 

X L 

= - 184258.5823 

km 

y'i 

= 329791.1351 

km 

z h 

= 103840.2349 

km 

Satellite position 

and velocity, in frame 

2 2 , at time t 

x i 

= - 4009.582237 

km 

x 2 

= + 103.9008135 

km 

X 3 

= - 5269.570696 

km 

*1 

= - 5.978210289 

km sec' 1 

x 2 

= + 1.710442216 

km sec" 1 

X 3 

= + 4.584971066 

km sec" ! 


= 398601 

km 3 sec' 

R 

= 6378.145 

km 


e 2 = 6.693421623E - 03 


k 20 - u.j 

k 2 , = 0.01 
k 30 = 0.1 
k 3] = 0.01 
k 22 = 0.1 

m' L = 7.3693281 E + 22 kg 

M = 5.97 31613 E + 24 kg 

a' L = 384400 km 
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Results 


0 = 0.1516444800 

deg 

£ = 0.1744119454 

deg 

z = 0.1744282487 

deg 

: 0 = 0.5000000017 

centuries 

T = 0.2723967010 

centuries 


X x = + 0.9999779632 
Y x = - 0.0060883617 
Z x = - 0.0026466801 
X y = +0.0060883617 
Y y = +0.9999814657 
Z y = - 0.0000080574 
X ? = +0.0026466801 
Y y = - 0.0000080567 
Z z = +0.9999964975 


\[ = - 186537.2414 

km 

y' L = + 328662.3531 

km 

z' L = + 103349.5407 

km 

r* = 391785.9267 

km 
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x* = - 188928.9046 


km 


y* = + 327293.3757 km 

z* = + 103349.5407 km 

X* = - 0.4822248369 
A<* = + 0.8353882910 
v* = +0.2637908451 
X = - 0.6054593779 
H = +0.0156893457 
v = - 0.7957215507 

r = 6622.380268 km 

a = 6567.451 187 km 

K. = 0.323073882E - 05 km 2 sec ’ 2 

n = 0.118624368E- 02 sec -1 

a = 0.1659246878E - 01 
a = 0.971 1750904 
= 0.9811480551 

A„ = +0.1978107925 
A' = - 0.3489996026 
\' 2 = - 0.1208534947E + 01 
A 3 = - 0.3816194918 
A' = +0.661 1033498 
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B', = - 0.1 179169837 
B' 2 = +0.2042749770 
B' 3 = +0.1748980753 
B; = +0.56091 18737 
Bj = - 0.0001311647178 
B; = - 0.4603145005 
B; = - 0.1594002275E + 01 

Pq = +0.06443748370 
P n = - 0.5995183587 
P', = - 0.09451271981 
P, = +0.3663349027 
P ' 2 = - 0.3272838250 
P 2 = - 0.009499261487 
P ' 3 = - 0.1190554808 
P 3 = + 0.4817770751 
P; = +0.2062472668 
P 4 = - 0.01248435049 

S' = +0.0005712740433 
S, = +0.131 1342628B + 01 
S', = - 0.0009896538092 
S 2 = - 0.03398098796 
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s' 3 = +0.001471593023 

5 3 = +0.1319815046 

S; = +0.0009131459348 

5 4 = - 0.2215028279 

Sj = - 0.2135318121E - 06 

5 5 = + 0.0172503888 
Sg = - 0.004239372772 

5 6 = - 0.2915005769 
S' 7 = - 0.01468033233 

5 7 = + 0.007558767081 

T' = +0.02363141132 
T, = - 0.2300019019 
T' 2 = - 0.04093817407 
T, = +0.005960068473 
T'j = - 0.002454126571 
T 3 = - 0.1963411384b + 01 
T; = + 0.9131459348b - 04 
T 4 = + 0.1762545737 
T' 5 = - 0.21353 18121E - 07 
T 5 = - 0.01372650615 
= + 0.02595922645 
T 6 = - O.I257338I35H + 01 
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T ’ 7 = +0.08989303178 
T ? = +0.03260345555 

V 0 = + 1.0000000000 
V, = +0.1742073243 
V, = - 0.1494101174 
V 3 = + 0.1899534786E- 02 
V 4 = - 0.03055341082 

C„ = - 0.0001622651725 
C, = + 0.1353296915E + 01 
C 2 = + 0.863 1523953E + 06 
C 3 = + 0.5505311134E+ 10 
C 4 = + 0.351 1367268E+ 14 

F = - 0.007748690189 
H = +0.1048876747 

p, , = + 0.2040473678 
p , 2 = - 0.004135328808 
p )3 = + 0.1119464221 
p 2| = +0.03700735158 
p 22 = - 0.004983233113 


km 3 sec -2 
km 4 sec -2 
km 5 sec -2 
km 6 secT 2 
km 7 sec ' 2 
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p 23 = + 0.1740434693 
p 31 = +0.3829649087 
p 32 = - 0.01385120787 
p 33 = + 0.2036553308 


- 0.1929747594E - 

09 

km 

sec 2 

+ 0.9604765343E - 

10 

km 

~ 2 
sec 

- 0.1824727285E - 

09 

km 

sec" 2 


T = 0.2824193799E - 09 km sec" 2 


DEVELOPMENT STATUS OF THE COMPUTER ROUTINES 


By early summer 1978, all five tide routines had been coded and made part, as 
operational subroutines, of the TERRA and CELEST computer program systems. Program 
checkout was accomplished (separately for each tide routine) by first selecting typical input 
data that would exercise the algorithm to a sufficient degree. Subsequently, the programmer 
would run the algorithm under test on the computer, while the author would perform 
identical calculations on a programmable electronic calculator. Agreement to ten significant 
digits of the computer results with those of the manual computation was taken as evidence 
for successful computer implementation of the algorithm. 

As indicated by the last paragraph of the INTRODUCTION, the two ocean tide 
algorithms may be expected to produce nearly identical perturbing accelerations because the 
tidal ocean surface of the second algorithm is a functionalization of the discrete 
representation of the ocean surface on which the first algorithm is based. Like the routines 
for the two air tides and that for the solid earth tide, each of the two ocean tides had 
been checked out individually. In the process, to facilitate the manual calculations, a simple 
ocean surface had to be selected for each of the two algorithms. Because of the different 
mathematical structure of the two algorithms, simplicity meant that each algorithm required 
its own individual surface. Thus, it was impossible to directly compare the resulting 
perturbing accelerations. This in turn meant that, during this phase of the checkout, the 
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mutual consistency of the two ocean tide routines could not be verified. Because the two 
routines are so complex, tins remained a source of concern despite the caution practiced 
while deriving them. The question was finally settled by running each algorithm, using the 
complete input data base. The resulting perturbing accelerations agreed with each other to 
20 percent. Considering the differences in the physical definitions of the tidal mass layers 
(one being a discrete point set, and the other being a continuous surface), this was 
interpreted as satisfactory evidence for the validity of the algorithms. 

For operational use. that ocean tide routine was selected which employs the discrete 
tide surface representation. Realizing that the second algorithm contains a more accurate 
physical model (by virtue of its continuous tidal mass layer as compared with the point 
mass approximation inherent in the first ocean tide routine), it was yet considered a 
decisive advantage of the algorithm chosen that it is capable of accepting among its input, 
and of being capable to directly process, without the need for any intervening 
computational steps, the earlier mentioned tide amplitude and phase tables that are in-house 
products of NSWC. This aspect is particularly important because it is expected that, during 
the useful life ot TFRRA/CLLhST, the present M, tide tables will be augmented by 
similar tables for various other constituents of the ocean tide spectrum. Once available, a 
mere extension of the present algorithm will make possible the immediate use of the 
additional program input, without first having to obtain the additional functionalizations of 
the tidal ocean that would be needed if the second ocean tide routine had been adopted 
for inclusion in TERRA/CELEST. 
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NOTES CONCERNING THE SECOND OCEAN TIDE ROUTINE 


As mentioned in the INTRODUCTION, the second computer routine for the ocean 
tide is based on a functionalization of the tidal ocean surface. This functionalization 
represents the height of the tidal surface above mean sea level, h, as an expansion in 
surface harmonics. 


h = f cos (ot* + x - 5) = ? cos 6 cos (ot* + x) + f sin 6 sin (at* + x) (A-l) 


15 n 

f cos 6 = £ £ ( C nm cosmX + S nm sin mX)P™ (sin 0) 

n = 0 m =0 


(A-2) 


1 5 n 

f sin 5 = V* V 1 (C' cos mX + S' sin mX)P m (sin 0) 

a—/ Z—< nm nm n v ' 

n = 0 m = 0 


(A-3) 


where 

X = longitude 
0 = latitude 

C„„, S„ m , C' , S' = expansion coefficients 

nm nm nm nm r 

a = rate of the mean mean longitude of the moon 
X = mean mean longitude of the moon, at the beginning 
of the day UT 

Note that h is essentially a function of longitude, latitude, and time. This function is valid 
everywhere on the globe (sea and land - see the comments at the end of this appendix). 

The surface density corresponding to the tidal mass layer is 


1 ( | 
s = ph = — cos (ot* + x) + (3 sin (at* + x) j 


(A-4) 


a = Gp f cos 5 


(A-5) 


(3 = Gp f sin 8 


(A-6) 
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The gravitational potential associated with the tidal bulge is then 


sdS' 

0 = Gff - = 0. cos (at* + x) + 0, sin (at* + x) 

|r- r | 1 1 


(A-7) 


</>, = 


. a dS' 

' It- r'l 


(A-8) 


(3 dS' 

>2 77 

z |r - r | 


ffdS 

I 


2 7r 

■' = R 2 j dX' f 

E 


+ 7T /2 

cos O' dd' 

0 *'-ir/2 

R n 


|F-7'| ~ 


P n (sin 0) P (sin 0') 


(A-9) 


(A-10) 


" (n - m)! , i 

+ 2£ ;--P™(sin0) P" 1 (sin0) cosm(X-X')f 

n “, (n + m)! 11 n > 


(A-ll) 


Now, evaluate 0. 


= GpR z fdX' I"cos0'd0' 

*0 -ir/2 (If-r 


(A-12) 


Let 


sin 6 = y (A-13) 

sin O' = y' (A -14) 
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ami consider the term ij from the expansion for f cos 6 , namely 


(f cos 5) u - (Cjj cos jX' + S 4j sin jX')PJ(y') (A-15) 

Then, 

(0i)jj = GpR 2 C ( . J dX' ldy' P^(y') cosjX' 

o -1 

•)£;£! [ p „<y> P.(y') 

( n=0 r L 

„ ^ (n - m)! m -i i 

+ 2 £, (7T^! P n (y) p n(y ) cos m(X - X')J j 

.+1 

dy' P J j(y') sin jX' 

i 

( 00 R n r 

* E ^ k(y) p n (y') 

I n = 0 1 L 

+ 2 n ^, (n + m)! C (y) P n <y'> c«sm(X- X')] J (A-16) 


+ GpR 2 S 


« I dX ' j 


Repeated application of the orthogonality relationships for the Legendre functions, as 
well as for the sine and cosine functions, makes possible the following equations: 


( 0 ,) 


>.j=o 


■ C„RiC, 0 2„ L' P,(y') f \ -51 

I ", F 

■ 2 ’ r2g “ 2m c i. p ,<y> 


P n (y) P n (/) 


(A-17) 
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- 2GpR I c j j t ± £ jgy> T(y') d/ 

* 1^' (y) fcosjX' los m(X - X') dX'j 

+ 2c ' R2s “ #o i, 

J 2n 

sin jX' cosm(X-X') dX'j 

I cos jX cos m(X - X ) dX r = 7r 6 cos iX 
J o 

C” . 

JsinjX cosm(X-X') dX' = 7r 5 mj sinjX 


• n = m 

if 

0 n # m 


(0 l >ij# 0 2jrGpR 2 prj pi( sin 0) * (C fj cos jX + S- sin jX) (A- 


and consequently, 


2 Rl 

01 )‘j = 2 " r2g P £7^ prn * (Cy cosjX + Sy sinjX) P](sin 


Thus, the following results in a similar fashion: 


2 R 1 

(0 2 >u = 27rR2Gp 2iTT * < c ij cosjX + S'j sinjX) Pj(sin 
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Noting that the original expansion of the tidal surface is limited to n x = 15, the 
tide potential is now 


0 = V V [F U + H V 1 

v\x vn v n vnj 

l> = 0 m=0 


(A-25) 


27rR 2 Gp 2 r . I 

F = - - C cos(at* + x) + C sin (at* + x) (A-26) 

vm 2 v + 1 *- "m J 

2rrR 2 Gp 2 r , 

H =-IS cos (at* + x) + S sin (at* + x) (A-27) 

I'm p, 2v + 1 L " M l ' M J v ' 

= r TTT p £ cos pX (A-28) 

= ( sin0 > sinpX (A-29) 


Up to here, each of the physical quantities occurring in the equations may still be 
entered in terms of arbitrary units. However, for the purpose of the second ocean tide 
computer algorithm, it is necessary to employ compatible units: kilometers for length, 
kilograms for mass, and mean solar seconds for time. As obtained from NOAA, the 
coefficients C^, S^, C' , and are listed in meters. To accommodate the existing 
computer card deck for these data, it is necessary to multiply the right sides of Equations 

A-26 and A-27 by 10“ 3 prior to their use in the main part of this report. 

It should finally be mentioned that the expansion for the tidal surface (Equations A-l 

through A-3) is intended to be valid for both sea and land in the sense that it is a fit to 

the NSWC/Schwiderski tidal heights and phases on the ocean and to zero tidal heights on 
land. Because of the relatively limited number of terms employed in the expansion, this fit 
is somewhat imperfect, especially near the coast lines. This imperfection is intentional. The 
reason for it is that the expansion was originally meant to be the basis for finding a tide 
potential for satellite work rather than a continuous, highly accurate description of the 
tidal surface. Satellite orbits appear to be much more sensitive to the low-order harmonics 
in the tidal disturbing potential than to the higher ones. Based on information from an 
outside source, we in fact anticipate that certain of the second- and fourth-order harmonics 
will be the main contributors to the orbital perturbation. Thus, the emphasis on the 
low-order terms apparent in the expansion for the tidal surface is rather an advantage than 
a defect as far as our own satellite work is concerned. 
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SOLID EARTH TIDE WITH THAT CAUSED BY THE OCEAN TIDE 
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A COMPARISON OF THE PERTURBING ACCELERATION DUE TO THE 
SOLID EARTH TIDE WITH THAT CAUSED BY THE OCEAN TIDE 


Preliminary runs of the first computer routine for the ocean tide and of the routine 
tor the solid earth tide, all of which used complete sets of input data (as opposed to the 
restricted input data of the numerical trial cases listed in the main body of this report), 
have resulted in the following perturbing accelerations. 

The satellite coordinates and the perturbing acceleration components listed are referred 
to the inertial frame corresponding to the mean equinox and mean equator of J = 1977. 
n = 202, and t* = 0 sec. 


Time t 


J = 1977 
n = 202 
t* = 50000 sec 

Satellite coordinates at t 

x, = 3151.52923 km 
x 2 = 5458.60875 
x 3 = 3639.07250 


Perturbing acceleration due to ocean tide at t 


T x - + 0.594629E - 11 km sec 2 
T y = + 0.134622E - 10 
T z = - 0.258193E- 11 


Perturbing acceleration due to solid earth tide at t 


T x = + 0.149368E - 09 km sec~ 2 
T y = + 0.125475E - 10 
T z = - 0.184573E - 10 

T x = - 0.493523E- 10 
T y = +0.381372E- 10 
T z = + 0.713397E - 12 


due to moon 


due to sun ' 


Nominal values for 
Love numbers 
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T x = + 0.164304K - 09 
= + 0.013802E - 09 
T z = -0.020303E - 09 

T v = -0.054288E - 09 
T y = + 0.041951 E- 09 
T z = + 0.000785 E- 09 

T x = + 0.17924IE - 09 
T = + 0.150570E - 10 
T z = - 0.221487E - 10 

T x = -0.592228E- 10 
T v = +0.457647E- 10 
T z = +0.856076E- 12 


The magnitudes of the perturbing accelerations are then 


TfOcean tide) - T 0( . EAN = 1.49E- 11 km sec 2 

T(Solid earth tide due to moon, for nominal Love Numbers) 

= T SOL(d = 1.51026E- 10 km sec' 2 

T(Solid earth tide due to moon, for Love numbers increased by 10 percent) 
= T SOL | L) + 10 = 1.66048E- 10 km sec' 2 

T(Solid earth tide due to moon, for Love numbers increased by 20 percent) 
= T so l i d + 2 o = 1«I-31E 10 km sec' 2 

Note that 



^SOLID+IO OLID ^ 'oCtAN (B-l) 

or 

^SOL1D + 10 ~ ^SOLID I ^ I^OCEAnI (B-2) 
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This indicates that, for the particular combination of time and satellite position 
chosen, if one increases the Love numbers by ten percent, the perturbing acceleration due 
to the lunar component of the solid earth tide increases by an amount roughly equal to 
the perturbing acceleration due to the ocean tide For the particular example chosen, the 
two equations are accurate to better than two percent, which is suspected to be a 
coincidence. But we are certain that this example reveals the order of magnitude of a 
relationship existing between the gravitational effects of the ocean tide and the solid earth 
tide. 
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